/******************************************************************************* 
                                SOLVE

  Title:        Single Particle solver.
  Author:       Hanna Szoke
  Date:         10/1/94
  Function:     This program estimates electron densities using A. Szoke's 
		holographic restoration method and Dennis Goodman's 
		non-linear complex conjugate gradient optimizer, CCG.

		Output includes a log file and optionally a listing of the
		cost function, plus the solution file in physical space 
		after each outer iteration, written as a binary file in units 
		of el/voxel.  This serves as input to the postprocessor, REGRID.
		Other optional output includes an outlier report, listing all
		reflections that deviate by more than 4*sigma from the fobs.

		Solve.c works in conjunction with the following source files:

		ccgutil.c is the interface to ccg.c;
		cost.c has the cost function calculations;
		fftshare.c is the interface to the fft package;
		highres.c for ungridded high-resolution points;
		hklutil.c, hklread.c and hklwrite.c handle (hkl) space work;
		init.c and basics.c handle initialization procedures;
		qshare.c handles general functions related to Np space;      
		crysutil.c  handles crystal symmetry;
		util.c has low-level worker routines.

*******************************************************************************/

#include "speden.h"		/* util.h, cellparams.h, symmetry.h, dims.h */
#include "speden_Np.h"
#include "cost.h"
#include "sp.h"

static	COMPLEX	*R0 ;		/* Model structure factors for known atoms */
static	real	AveWt ;

real	sp_min ;	/* for comparing cost */
real	hkl_cost ;	/* current value of hkl cost function */

FILE	*fp_cost ;					/* in util.h */
char	cost_filename[MAXSTRING] ;			/* in cost.h */
char	costfile_option = FALSE ;			/* in util.h */
char	*maskfc ;	/* mask in speden.h identifying fc input values */

COMPLEX *R ;    	/* Model structure factors for known+found atoms */
real	norm_fac ;	/* 1/sqrt(AveWt), in cost.h */
int	Nfft ;		/* Nfob or (Nfob+Nfc) if there is phase extension */

void	convolve(int, real *, real *, COMPLEX *) ; 
static	void	solver_setup(char *) ;
static	real	getRfactor() ;
static	void    makeStructureFactor() ; 
static	void	outliers_summary() ;
static	void	outliers_details() ;
static	void	init_outliers_details() ;
static	void	calcChiSq() ;
static	real	getDiscrep() ;
static	void	rereadR0() ;
static	void	report_costs(char *) ;
static	void	report_mode() ;
static	void	report_memory() ;
static	int	setup_fobs() ;
static	void	setup_fcalc() ;
static	void	alloc_hklArrays() ;

void	solve_main(int argc, char *argv[])
{
     int	checkSolutions() ;		/* in 'qshare.c' */
     void	writeSol(char *) ;	
     void       updateKnownValues() ; 
     void       reinitNump() ;
     int       	HoloMethod() ;			 /* in 'ccgutil.c' */
     real	doStdDev(char *) ;
     void	reportTotFunctCalls() ;
     void	do_inplace_cs_symmetrization(real *, real *) ;
		      			/* in 'crysutil.c' */

     void	hrWriteLists(char *, real *) ;
     void	hrSymmetrize(real *) ;

     real	R_fac ;			/* Current R factor */
     int	n ;
     real	R_fac_init ;		/* Initial R factor */
     int        iter ;
     real	stdev_init, stdev_uns, stdev_end ;	/* ret'd by doStdDev */
     int	stdev_crit = 0;     /* continue? based on standard deviation */
     int	Holo_crit ;		/* continue? ret'd by HoloMethod */
     int	update_crit ;		/* continue? ret'd by checkSolutions */
     char	problemName[MAXSTRING] = "" ;
     char	out_name[MAXSTRING] ;	/* generic name, used for output */

#define DEL_STD_DEV	0.005	/* min. fract. change in std. dev. over an outer 
				iteration before code quits */

  /******************************************************
  Check for right ballpark (non-fatal).
  ******************************************************/

     if (argc > optind+2) { 	/* too many arguments */
          sprintf(message, "\n\n\tIgnoring argument(s) %s ...", argv[optind+2]) ;
	  printTwice(message) ;
     }

  /******************************************************
  Pull off any prefix defining whereabouts of input file
  and remove the .inp extension.
  Fetch input, prepare for doing optimization.          
  ******************************************************/

     strip_path(input_filename, pwd_name) ;
     strncpy(problemName, pwd_name, strlen(pwd_name) - 4) ;
     sprintf(out_name, "%s", problemName) ;

     solver_setup(problemName) ;

     makeStructureFactor() ;

  /******************************************************
  Free R0 to save space; it will be re-read, used & freed
  within the main loop after completing HoloMethod().
  Prepare an initial outlier report.  
  ******************************************************/ 

     e_free(Nfft*sizeof(COMPLEX), R0) ; 

     for (n = 0; n < Nfob; n++) 
	  *(ReList + n) = Cabs(R+n) ;
  
     if ((R_fac_init = R_fac = getRfactor()) < R_stop) {
          sprintf(message, 
	       "\nQuitting immediately - initial R factor is less than %g",
	       R_stop) ;
          Warning(message) ;
          return ;
     }
        
     outliers_summary() ; 
     if ((verbose) && (strcmp (md_filename, "empty") != 0)) 
           init_outliers_details() ; 

     if (useSig)
	  calcChiSq() ;

  /************************************************************
  Iterate; 50 is a magic number to stop code if it doesn't 
  converge, but we've never hit it.  Typically, iter reaches
  2 to 5 before the code quits.
  ************************************************************/

     for (iter = 0;iter < 50;iter++) { 

          printTwice("\n") ;
          printTwice(timestamp()) ;
          printTwice("Applying holographic reconstruction - ") ;

	  sprintf(message, "iteration # %d", iter+1) ; 
	  printTwice(message) ;

	  if (costfile_option)
	       fprintf(fp_cost, "\niteration # %d\n", iter+1) ; 

       /************************************************************
       Do holographic reconstruction, symmetrize.  Add new 
       contributions to snump, write out an electron density map.
       ************************************************************/

          sp_min = getDiscrep() ;

	  stdev_init = doStdDev("Initial standard deviation") ;  
          if (iter == 0)
	       report_costs("Initial") ;
						/**************/
          Holo_crit = HoloMethod();		/* THIS IS IT */
						/**************/
	  update_crit = checkSolutions(); 
		
          if (Ncs > 1) {

             /** below is MISUSE of snump! If we ever start using non-P1
                 crystals, do_inplace_cs_symmetrization() must be
                 revisited wrt use of nump (csol) and snump.   **/

	     stdev_uns = doStdDev("Standard deviation before symmetrization") ;

             do_inplace_cs_symmetrization(nump, snump) ;  
	     if (HighRes)
	          hrSymmetrize(nump+Nptotal) ;

	     stdev_end = doStdDev("Standard deviation after symmetrization") ; 
          }
	  else
	     stdev_end = doStdDev("Standard deviation") ; 

       /************************************************************
       If over one outer iteration, the standard deviation does not 
       decrease by more than DEL_STD_DEV (3%), quit.     
       *************************************************************/

          if ((stdev_init - stdev_end) < DEL_STD_DEV * stdev_end) {
	     printTwice("\nStopping - standard deviation is not decreasing ...\n") ;
             stdev_crit = STOP ;
          }

          updateKnownValues() ; 

          writeSol(out_name) ; 
          
          if (HighRes) {
               sprintf(message, "%s.list", problemName) ;
               hrWriteLists(message, totp + Nptotal) ;
	       if (verbose) {
                    sprintf(message, "it%1d.list", iter+1) ;
                    hrWriteLists(message, totp + Nptotal) ;
	       }
          }
       /************************************************************
       Do an analysis of the current solution in (hkl) space.
       ************************************************************/

	  rereadR0() ;
          makeStructureFactor() ;
          e_free(Nfft*sizeof(COMPLEX), R0) ; 

          for (n = 0; n < Nfob; n++) 
	       *(ReList + n) = Cabs(R+n) ;

          R_fac = getRfactor() ;

          outliers_summary() ; 
          if (verbose) 
               outliers_details() ;
          
	  if (useSig)
	       calcChiSq() ;

       /************************************************************
       Decide whether to go on, based on various figures of merit.
       ************************************************************/

          if (R_fac < R_stop) {
	       sprintf(message, 
		 "Stopping - Rfac is less than %g\n", R_stop) ;
               printTwice(message) ;
	       break ;
          }

          if ((stdev_crit == STOP) || 
	      (Holo_crit  == STOP) ||
	      (update_crit == STOP)) 
	       break ;

          if (oneIt)
               break ;

          reinitNump() ;
	  report_costs("Current") ;
     }

    report_costs("Final") ;

    sprintf(message, 
	 "\nOverall R factor changed from %g to %g", R_fac_init, R_fac) ;
    printTwice(message) ;

    reportTotFunctCalls() ;
}

void	solver_setup(char *problemName)
{
     int        HDFdatalength(char *, char *) ;
     int        read_HDFfcalc(char *) ; 
     int        readFcList(char *) ;    
     void	allocNpArrays() ;		/* in 'qshare.c' */
     void	alloc_fc(int);			/* in 'hklread.c' */
     void	initNparrays() ;		/* in 'qshare.c' */
     void       echo() ;			/* in 'init.c' */
     void       echo_ignored() ;		/* in 'init.c' */
     void       readInput() ;			/* in 'init.c' */
     void	initfft(int, int) ; 		/* in 'fftshare.c' */
     int	k ;

	/* set up cost file */

     hello(caller) ;

     if (verbose)
	  costfile_option = TRUE ;

     if (costfile_option) {
	  sprintf(cost_filename, "%s.cost", problemName) ;

          if ((fp_cost = fopen(cost_filename,"w")) == NULL) 
          {
               sprintf(message, "Cannot open %s", cost_filename);
               Fatal(message) ;
          }
     }

  /************************************************************
   Fetch input conditions, report.                         
  ************************************************************/
 
     readInput() ;

     report_mode() ;
     echo() ;
     echo_ignored() ;

  /************************************************************
   Do fobs part of the hkl-space setups.  Then allocate and 
   initialize arrays for real-space arrays.  Finally, generate
   fcalc files.  
  ************************************************************/

     Nfob = setup_fobs() ;
     report_memory() ;

     allocNpArrays() ;
     initNparrays() ;
     
     sprintf(message,
           "\n\t\tThis problem has %d equations, %d unknowns\n", Nfob, 
		Npextended) ;
     printTwice(message) ;

     for (Nfc = 0, k = 0; k < Nconstraints; k++) 
	
          switch (con_type[k]) {

	  case PHASE_EXT :
              if (HDFIN) {
                    Nfc = HDFdatalength(fc_filename[k], "/fcalc/Nfc") ;
                    alloc_fc(Nfc) ;
                    printTwice("Reading an HDF file for phase extension") ;
                    Nfc = read_HDFfcalc(fc_filename[k]) ;
              }
              else {
                    Nfc = filelength(fc_filename[k]) ;
                    alloc_fc(Nfc) ;
                    printTwice("Reading a list file for phase extension") ;
                    Nfc = readFcList(fc_filename[k]) ;
              }
              Nfft = Nfob+Nfc ;
              initfft(3, Nfft) ;
              break ;

	  default :
	       break ;
          }

     if (Nfc == 0){
         Nfft = Nfob ;
         initfft(2, Nfft) ;
     }

     alloc_hklArrays() ;
     setup_fcalc() ;
}

int	setup_fobs()
{
	/* functions in 'hklutil.c' or hklread.c */

     int	Nlen ;
     int        HDFdatalength(char *, char *) ;
     void	SigmaRange(real *, real *, real *, int) ;
     int	checkSigmas(real *, int) ;
     real	getNormFac() ;
     void	fsetDsqLimit(int, char) ;
     int	readfobs(char *) ;
     int	read_HDFfobs(char *) ;
     void	scalefobs(real *, real) ;

     /* local variables */

     int	n ;
     real	sigmin, sigmax ;	/* unused here */
     real	*wt ;
     real	*sarray ;

  /************************************************************
   Set up Fo arrays in (h,k,l) --  reciprocal lattice triplets   
  ************************************************************/
 
     if (HDFIN)
          Nlen = HDFdatalength(fo_filename, "/fobs/Nfob") ;
     else
          Nlen = filelength(fo_filename) ;

     printTwice("Setting up initial arrays ...") ;

     sprintf(message, "setup_fobs") ;

     Folist =  (real *) e_malloc(Nlen*sizeof(real), message) ;
     Siglist = (real *) e_malloc(Nlen*sizeof(real), message) ;
     Hlist =  (float *) e_malloc(Nlen*sizeof(float), message) ;
     Klist =  (float *) e_malloc(Nlen*sizeof(float), message) ;
     Llist =  (float *) e_malloc(Nlen*sizeof(float), message) ;
     hkl_weight = (real *) e_malloc(Nlen*sizeof(real), message) ;  

  /************************************************************
  Read diffraction pattern, Folist, for the full molecule.  
  Note that Nfob is the actual number of fobs used (i.e., it
  excludes values at a resolution too high for the current 
  problem) and it includes the (000) term from the input file
  as its first entry.
  ************************************************************/

     printTwice("\nReading fobs file ") ;

     if (HDFIN)
          Nfob = read_HDFfobs(fo_filename) ;
     else
          Nfob = readfobs(fo_filename) ;

     scalefobs(Folist, fscale) ;
     scalefobs(Siglist, fscale) ;

     fsetDsqLimit(Nfob, TRUE) ;

  /************************************************************
  Check sigma range and reset flag if there were missing sigmas 
  in input.  Set criterion for discrepancy-criterion-satisfied 
  and normalize the weights.
  ************************************************************/

     useSig = checkSigmas(Siglist, Nfob);

     if (useSig) {
          printTwice("\nFor Fobs list:") ;
          SigmaRange(&sigmin, &sigmax, Siglist, Nfob) ; 
     }

  /************************************************************
  Handle weights, including normalization.
  ************************************************************/

     for (n = 0, wt = hkl_weight; n < Nfob; n++, wt++) 
          *wt = 1 ;
     
     *hkl_weight = 1./ SQRT2  ;
     
     AveWt = getNormFac() ;
     norm_fac = 1. / sqrt(AveWt) ;

     for (n = 0, sarray = Siglist, wt=hkl_weight; n < Nfob; n++, sarray++, wt++) 
	 if (*sarray > 0)
	       *wt *= norm_fac / *sarray ;
         else
	       *wt = 0 ;

     fprintf(fp_log, "Hkl weight normalization factor = %g\n", norm_fac) ;

     return (Nfob) ;
}

void	setup_fcalc()
{
     void	describeFR(char *) ;			/* in hklutil.c */
     char	*get_P1mask() ;				/* in hklutil.c */
     void	hrPrepareSFT() ;			/* in highres.c */
     void	hrDoSFT(real *,COMPLEX *) ;		/* in highres.c */
     int	n ;
     real       delta, phdr ;

  /************************************************************
   Set up arrays in (h,k,l) --  reciprocal lattice triplets   
  ************************************************************/
     expfac = (real *) e_malloc(Nfft*sizeof(real), "setup_fcalc") ;

     for (n = 0; n < Nfob; n++){
          delta = PISQ * eta * dr * dr ;
          *(expfac+n) = exp(-(delta * fo_nusq(n))) ;
     }
     if (thereis_phext) {
         phdr = dr * phase_ext_res / input_res    ;
         delta = PISQ * eta * phdr*phdr ;
         for (n = 0; n < Nfc; n++)
              *(expfac+Nfob+n) = exp(-(delta * fc_nusq(n))) ;
     }

     if (HighRes)
          hrPrepareSFT() ;

  /************************************************************
  Generate a diffraction pattern, R0, corresponding to the known 
  part of the molecule, from knownp; this will be 
  added into an iterated R within the loop. Enforce consistent 
  masking.  Provide some diagnostics on Folist and R0. 
  ************************************************************/

     printTwice("Generating starting position in reciprocal space.") ;

     maskfc = get_P1mask() ;
     convolve(Nfft, knownp, expfac, R0) ;

     if (HighRes)
          hrDoSFT(knownp + Nptotal, R0) ;

     describeFR(maskfc) ;
}

void alloc_hklArrays()
{
  /************************************************************
   Arrays in hkl (Native); Oar is local to cost.c.
   Note that R0 is allocated, used and then freed;
   it will be brought in again and used in rereadR0().
  ************************************************************/

     strcpy(message, "alloc_hklArrays" );

     maskfc = (char *) e_malloc(Nfob*sizeof(char), message) ;
     
     R =  (COMPLEX *) e_malloc(Nfft*sizeof(COMPLEX), message) ;
     R0 = (COMPLEX *) e_malloc(Nfft*sizeof(COMPLEX), message) ;
     ReList = (real *) e_malloc(Nfob*sizeof(real), message) ; 
     Oar = (COMPLEX *) e_malloc(Nfft*sizeof(COMPLEX), message) ; 
     Phases = (real *) e_malloc(Nfft*sizeof(real), message) ; 
}

				/****************************************
                 		 Calculate the R factor for each part of
				 the problem (native plus derivatives) 
				****************************************/

real getRfactor()
{
     void	shellHeader() ;	/* in hklutil.c */
     real	get1Rfac(real *, real *, int) ;		

     real	thisr ;

     shellHeader() ;
                                                                     
     thisr = get1Rfac(Folist, ReList, Nfob) ;

     return(thisr) ;
}

     		/************************************************************
     		Assemble structure factors, R, from FFT-generated reflections 
     		from found atoms plus R0 (reflections from known atoms).
     		************************************************************/

void makeStructureFactor() 
{
     void	hrDoSFT(real *, COMPLEX *) ;

     static int first_time = TRUE ;
     int        n ;

     /************************************************************
     Make structure factors for the grid(s).  
     Propagate R throughout array.
     ************************************************************/

     convolve(Nfft, snump, expfac, R) ;

     if (first_time)
          first_time = FALSE;
     else
          if (HighRes) 
               hrDoSFT(snump+Nptotal, R) ;

     /************************************************************
     Add in known part of the native structure.
     ************************************************************/

     for (n = 0; n < Nfft; n++) {
          (R+n)->re += (R0+n)->re ;
          (R+n)->im += (R0+n)->im ;
     }
     
}

void	rereadR0() 
{

     void	hrDoSFT(real *, COMPLEX *) ;			/* in highres.c */

  /************************************************************
  Regenerate  R0 in (h,k,l) --  reciprocal lattice triplets.
  Get the diffraction pattern, R0, corresponding to the known 
  part of the molecule.  
  ************************************************************/

     R0 = (COMPLEX *) e_malloc(Nfft*sizeof(COMPLEX), "rereadR0") ;

     convolve(Nfft, knownp, expfac, R0) ;
     if (HighRes)
         hrDoSFT(knownp + Nptotal, R0) ;
}

				/*********************************************
	        		Calculate summary information about the
				outliers for a data set. 
				*********************************************/
void outliers_summary()
{
     real	spread ;
     real	buckets[5] ;	/* we use 5 buckets for collecting stats.  */
     real	p1, p2, p3, p4 ;
     real	numer ;

     int	Ntot = 0 ;
     int	n, b ;
     char	gotit ;
     char	which[10] ;

     if (useSig)
          printTwice(
	  "\n\tPercentages of Fcalcs within sigma intervals of Fobs:\n") ;
     else
          printTwice(
	  "\n\tPercentages of Fcalcs within intervals of Fobs:\n") ;
     printTwice(
     "interval:                1          2          3          4   outliers") ;

     for (b = 0; b < 5; b++) 
       buckets[b] = 0 ;

     strcpy(which, "native") ;
	       
     for (n = 0; n < Nfob; n++) {

            numer = fabs(*(Folist+n) - *(ReList+n)) ;

            spread = (useSig) ? numer / *(Siglist+n) : numer ;
            Ntot++ ;

            for (b = 0; b < 4; b++) {
	         gotit = FALSE ;
	         if (spread < (real) (b+1)) {
		      buckets[b] += 1. ;
		      gotit = TRUE ;
		      break ;
                 }
            }
            if (!gotit) 
	         buckets[4] += 1. ;
     }

     for (b = 0; b < 5; b++) 
         buckets[b] *= 100. / (float) Ntot ;

     sprintf(message, "%% %6s:       %10g %10g %10g %10g %10g", 
     which, buckets[0], buckets[1], buckets[2], buckets[3], buckets[4]) ;
     printTwice(message) ;

     p1 = erfc(1./SQRT2) ;
     p2 = erfc(2./SQRT2) ;
     p3 = erfc(3./SQRT2) ;
     p4 = erfc(4./SQRT2) ;

     sprintf(message, "%% for Gaussian: %10g %10g %10g %10g %10g", 
     100*(1.-p1), 100*(p1-p2), 100*(p2-p3), 100*(p3-p4), 100*p4) ;
     printTwice(message) ;
}

				/*********************************************
	        		Print out details about all outliers in
				a data set at start of run.
				*********************************************/

void init_outliers_details()
{
     FILE	*fp ;
     char	filename[12] ;
     real	spread ;
     real	rrr;
     real	numer ;
     real	srat ;

     int	n ;

     sprintf(filename, "outlier00") ;
     if ((fp = fopen(filename, "w")) == NULL)
          Fatal("Cannot open outlier file") ;

     if (useSig) {
          fprintf(fp, "    Outliers whose weighted deviation > sigma\n\n") ;

          fprintf(fp, 
"    Res (nm) deviation      fh     fk     fl        Fobs     |Fcalc|         Sig     dev/Fobs \n") ;
     }

     else {
          fprintf(fp, "    Outliers whose weighted deviation > 1\n\n") ;

          fprintf(fp, 
"    Res (nm) deviation      fh     fk     fl        Fobs     |Fcalc|     dev/Fobs \n") ;
     }

     for (n = 0; n < Nfob; n++) {	

         numer = (*(Folist+n) - *(ReList+n)) ;
         spread = (useSig) ? numer / *(Siglist+n) : numer ;

	 if (fabs(spread) >= 1.) {
         
	    srat = spread / *(Folist+n) ;
   	    rrr = 1. / sqrt(fo_nusq(n)) ;

	    if (useSig)
	              fprintf(fp, 
	 " %10.4g %10.4g  %6.2f %6.2f %6.2f %11.4g %11.4g  %10.4g  %10.4g\n",
          rrr, spread, *(Hlist+n), *(Klist+n), *(Llist+n),
	  *(Folist+n), *(ReList+n), *(Siglist+n), srat) ;

	    else
	              fprintf(fp, 
	 " %10.4g %10.4g  %6.2f %6.2f %6.2f %11.4g %11.4g  %10.4g\n",
          rrr, spread, *(Hlist+n), *(Klist+n), *(Llist+n),
	  *(Folist+n), *(ReList+n), srat) ;

         }
    }
    fclose(fp) ;
} 

				/*********************************************
	        		Print out details about the far outliers in
				a data set. 
				*********************************************/
void outliers_details()
{
     FILE	*fp ;
     char	filename[12] ;
     real	spread ;
     real	rrr;
     real	numer ;

     int	n, b ;
     char	gotit ;

     sprintf(filename, "outlier0") ;
     if ((fp = fopen(filename, "w")) == NULL)
     Fatal("Cannot open outlier file") ;

     if (useSig) {
          fprintf(fp, 
  "    Outliers whose weighted deviation > 4*sigma\n\n") ;

          fprintf(fp, 
"   fh     fk     fl        Fobs     |Fcalc|         Sig     Deviation     Res (nm)\n") ;
     }

     else {
          fprintf(fp, 
  "    Outliers whose weighted deviation > 4 \n\n") ;

          fprintf(fp, 
"   fh     fk     fl        Fobs     |Fcalc|     Deviation     Res (nm)\n") ;
     }

     for (n = 0; n < Nfob; n++) {	

         numer = (*(Folist+n) - *(ReList+n)) ;

         spread = (useSig) ? numer / *(Siglist+n) : numer ;

         for (b = 0; b < 4; b++) {
	      gotit = FALSE ;
	      if (fabs(spread) < (real) (b+1)) {
	           gotit = TRUE ;
	           break ;
              }
         }
         if (!gotit) {

		 rrr = 1. / sqrt(fo_nusq(n)) ;
		 if (useSig)
	              fprintf(fp, 
	 " %6.2f %6.2f %6.2f %11.4g %11.4g  %10.4g %10.4g %10.4g \n",
	          *(Hlist+n), *(Klist+n), *(Llist+n),
		  *(Folist+n), *(ReList+n), *(Siglist+n),
                  spread, rrr) ;
		 else
	              fprintf(fp, 
	 " %6.2f %6.2f %6.2f %11.4g %11.4g %10.4g %10.4g \n",
	          rrr, spread, *(Hlist+n), *(Klist+n), *(Llist+n),
		  *(Folist+n), *(ReList+n)) ;

         }
    }
    fclose(fp) ;
} 

void report_mode()      /* Summarize Speden mode of operation.         */	
{
     int	k ;

     printTwice(
	"\t******************************************************\n") ;
     printTwice( "\tThis is a single-particle Speden run.") ;

     if (HighRes) {
          sprintf(message, "\tHigh Resolution processing is in effect.") ;
          printTwice(message) ;
     }
     if (Nconstraints == 0) {
	  sprintf(message, 
	       "\tThere are no Np constraints in the cost function.") ;
          printTwice(message) ;
     }
     else if (Nconstraints == 1) {
	  sprintf(message, 
	       "\tThere is 1 Np constraint in the cost function:") ;
          printTwice(message) ;
     }
     else {
	  sprintf(message, 
	       "\tThere are %d Np constraints in the cost function:",
	  Nconstraints) ;
          printTwice(message) ;
     }
     for (k = 0; k < Nconstraints; k++) 
	
          switch (con_type[k]) {

	  case TARGET :
	       sprintf(message, 
	            "\t   A stabilizing target with relative weight %g.", 
		    relwt_con[k]) ;
               printTwice(message) ;
	       break ;

	  case PHASE_EXT :
	       sprintf(message, 
	            "\t   A phase extension target at a resolution of %g.", 
		    phase_ext_res) ;
               printTwice(message) ;
	       sprintf(message, 
		    "\t   with relative weight %g", relwt_con[k]);
               printTwice(message) ;
               printTwice("\t   using list input in reciprocal space") ;
	       break ;

	  case CS :
	       sprintf(message, 
                    "\t   Crystal symmetry enforced with relative weight %g.", 
	            relwt_con[k]) ;
               printTwice(message) ;
	       break ;

	  default :
	       break ;

          }

     printTwice(
	  "\n\t******************************************************") ;
}

			/************************************************
       			Determine ChiSq 
			************************************************/
void	calcChiSq() 
{
     int	n ;
     int	Ntot ;
     real	chi ;
     real	chisq ;	
     real	wtchisq ;	

     Ntot = 0 ;
     chisq = 0 ;
     wtchisq = 0 ;

     for (n = 1; n < Nfob; n++) {	/* skip (000) term */

            chi = (*(Folist+n) - *(ReList+n)) / *(Siglist+n) ;
            chisq += chi*chi ;
            wtchisq += chi*chi / sqrt(1. + chi*chi/16.) ;
            Ntot++ ;
     }
     chisq /= Ntot ;
     wtchisq /= Ntot ;

     printTwice("") ;

     sprintf(message,
"Goodness-of-fit: chisq = %g, robust chisq = %g", chisq, wtchisq) ;
     printTwice(message) ;
}

  	/************************************************************
	Get "discrepancy value" (see paper VI, Appendix D).  
	************************************************************/

real getDiscrep() 	
{
     real	SsqWts ;	
     int	n ;
     real	*wt ;
     real	fedmin ;

     if (!useSig)
	  fedmin = DISCRP ;		/***	DISCRP defined in speden.h */

     else {
	  wt = hkl_weight ;

          for (n = 0, SsqWts = 0; n < Nfob; n++, wt++) 
	       SsqWts += *wt * *wt  ;
   
	  fedmin = 0.5 * SsqWts / AveWt ;

          sprintf(message, "%6.1e", fedmin) ;
          sscanf (message, "%g", &t1) ;
	  fedmin = t1 ;
     }
     fedmin *= discrp_frac ;            /* user-defined fraction */

     sprintf(message,
       "\nStopping criterion for the (hkl) cost function is %6.1e", fedmin) ;
     printTwice(message) ;

     return(fedmin) ;
}

void	report_memory()
{

     real	mem_Npspace, mem_recspace, mem_tot ;
     int	Nptot_gridPt ;

	/*******************************************************
	In Np space, there are 6 Speden real arrays (snump, 
	totp, knownp, nump, minp and maxp) plus 1 Speden int 
	array (typep) and 4 CCG real arrays (xn, gn, go and d),
	each with Nptotal elements.  For each target, there are 
	another 2 real arrays (target, weight) disregarding
	exotic constraints.  There are 2 complex FFT arrays in 
	Nptotal.  With phase extension, there is 1 more real array.
	*******************************************************/
     
     Nptot_gridPt = (10 + 2*Ntargets) * sizeof(real) + 1 * sizeof(int) ; 

     mem_Npspace = Nptot_gridPt*Nptotal + 2*sizeof(COMPLEX)*Nptotal ;

	/*******************************************************
	In reciprocal space, there are 3 Speden complex arrays 
	(R, R0 and Oar), 1 complex array (fftarray) plus 1 real 
	array (expfac), all in Nfft.  There are 4 real arrays in 
	Nfob (Folist, Siglist, ReList and hkl_weight) and 3 
	float arrays in Nfob (Hlist, Klist and Llist).  
	*******************************************************/

     mem_recspace = 5 * Nfft * sizeof(COMPLEX) +
		    1 * Nfft * sizeof(real) +
		    4 * Nfob * sizeof(real)  +
		    3 * Nfob * sizeof(float)  ;

     mem_tot = ((mem_Npspace + mem_recspace)*1e-9) + 0.5 ;

     if (mem_tot > 1.) 
          sprintf(message,
     "Approximate memory requirements in Gbytes: %6.2g", mem_tot) ;
     else {
          mem_tot = ((mem_Npspace + mem_recspace)*1e-6) + 0.5 ;
          sprintf(message,
     "Approximate memory requirements in Mbytes: %6.2g", mem_tot) ;
     }
     printTwice(message) ;
}
void report_costs(char *which)
{
        sprintf(message,
          "%s value of the (hkl) cost function is %6.1e", which, hkl_cost) ;
        printTwice(message) ;
}
